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Abstract 

The Alekseev-Grobner lemma is combined with the theory of modified 
equations to obtain an a priori estimate for the global error of numerical 
integrators. This estimate is correct up to a remainder term of order h 2p , 
where h denotes the step size and p the order of the method. It is applied 
to a class of nonautonomous linear oscillatory equations, which includes 
the Airy equation, thereby improving prior work which only gave the 
hP term. 

Next, nonlinear oscillators whose behaviour is described by the Emden- 
Fowler equation y" + t v y n = are considered, and global errors com- 
mitted by Runge-Kutta methods are calculated. Numerical experiments 
show that the resulting estimates are generally accurate. The main con- 
clusion is that we need to do a full calculation to obtain good estimates: 
the behaviour is different from the linear case, it is not sufficient to look 
only at the leading term, and merely considering the local error does not 
provide an accurate picture either. 

1 Introduction 

This paper is about the numerical solution of ordinary differential equations. 
Unfortunately, the outcome of the calculation is almost invariably not exact, 
but carries an error. The error committed by a single step of the numerical 
method is called the local (truncation) error. But in general, a single step is not 
sufficient to cover the time interval of interest, so we are interested in the error 
at a specific time, after several steps. This quantity is called the global error. 

It is important to have at least a rough idea of the size of the global error, 
so we need to estimate it. We can distinguish two approaches. One is to solve 
the equation numerically, and to use the result somehow to find an a posteriori 
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estimate. The other approach is to use the analytic solution, or at least some 
knowledge about it. This gives rise to a priori estimates. 

Both types of estimates have their respective strengths. If one has actually 
computed a numerical solution and wants to know how far off it is, one probably 
should use an a posteriori estimate. However, for the purpose of comparing 
different methods, or that of devising methods which are particularly suited for 
a certain class of problems, there is often no choice but to use a priori estimates. 

Skeel [p5| wrote a nice overview of various methods to obtain a posteriori 
estimates of the global error, but here we will concentrate on a priori estimates. 
It is rather hard to find accurate estimates of the global error, as the following 
quote by Lambert shows: [^0|, p. 57] 

The [local truncation error] and the starting errors accumulate to 
produce the [global error] , but this accumulation process is very com- 
plicated, and we cannot hope to obtain any usable general expression 
for the [global error]. However, some insight into the accumulation 
process can be gleaned by looking at an example. 

So the challenge is to find a class of equations which is specific enough to allow 
the resulting estimates to be of some use, yet general enough to encompass many 
problems arising in practice. The class we are aiming for here, can loosely be 
described as nonautonomous oscillatory equations. 

Among the first a priori estimates were those obtained by Henrici ]H| and 
Gragg ||, who expanded the global error in powers of the step size. This 
expansion was generalized and simplified by Hairer and Lubich Jl2] |, and used 
by Cano and Sanz-Serna || to analyse the behaviour of numerical methods 
for systems with periodic orbits. 

A parallel path of research is based on the Alekseev-Grobner lemma, to be 
introduced in the next section. We merely point to the work of Iserles and 



Soderlind 19 , Calvo and Hairer ||, and Viswanath pq . A comparison of both 
approaches can be found in the paper by Hairer and Lubich |b|| . 

The recent paper of Iserles gives an a priori estimate for the global error 
committed for certain linear oscillators. This estimate only concerns the first 
term in the expansion of the global error in powers of the step size. The goal 
of the research described here, is to extend this work in two directions. Firstly, 
more (but not all) terms of the expansion are derived. Numerical experiments 
will show that it is sometimes necessary to include these terms. Secondly, the 
effects of nonlinearity are studied by considering a specific example. 

This example is the Emden-Fowler equation y" + t v y n = 0. For a particular 
range of parameters, the solutions of this nonlinear equation are oscillatory. It 
was selected because it is neither linear nor autonomous, and the oscillatory 
behaviour of its solutions will probably lead to substantial cancellations causing 
more difficulties. Yet it is still amenable to analysis, in particular because we 
can derive its asymptotic solution for large t. Finally, it seems a logical sequel 
to the Airy equation y" + ty = 0, considered by Iserles Q. 

The remainder of this paper is organized as follows. The next section in- 
troduces the Alekseev-Grobner lemma and the theory of modified equations, 
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and shows how they can be used to estimate the global error. The resulting 
a priori estimates are exact up to a term of order h 2p , where h is the step size 
and p the order of the method. In Section ||, we show how these estimates 
look like for a class of linear oscillators, which includes the Airy equation. We 
condense the calculations because many are already described by Iserles Q. 
The Emden-Fowler equation is introduced in Section |4|, where we also describe 
its asymptotic solution. Section || contains the actual calculation of global error 
estimates. The estimates are supported by numerical experiments, as described 
in Section 0. The last section discusses the results, draws them into comparison 
with the work done by others, and presents some pointers for further research. 



2 Theory 

Suppose that we are solving the ordinary differential equation 

y'(t)=f(t,y(t)), y(*o)=yo, (1) 

where y(t) S R d is a vector. We assume that f is C°° to avoid technicalities 
which are irrelevant to our argument. 

Associated to this differential equation is the flow map, which is denoted 
here by $*(yo), where s stands for the initial time, t the final time, and yo the 
starting position. The map <&* : R d — > R d is defined by 

$*(y)=y and (y) = f (t, $*(y)) . (2) 

Its Jacobian matrix will be denoted by -0$* and is called the variational flow. 
It plays an important role in this paper, since it describes how perturbations are 
propagated by the flow of ([!]). This follows from the Alekseev-Grobner lemma 
(see eg. which extends the variation-of-constants technique to nonlinear 

equations. 

Lemma 1 (Alekseev-Grobner). Denote by y and z the solutions of 

y' = f(i,y), y(io) = yo> (3a) 

z' = f(t,z)+g(t,z), z(t ) = yo, (3b) 

respectively and suppose that di / dy exists and is continuous. Then the solutions 



of (3a) and of the "perturbed" equation (3b) are connected by 

<t)=y{t)+l D^ s (z(s))g(s,z(s))ds. (4) 

Jto 

Now suppose that one uses a numerical one-step method with fixed step 
size h to solve the differential equation (Q). Letting y = y(io)j we denote the 
subsequent results of the method by yi, ya, ya, . . . , and set tk = t + kh. 

Of course, we hope that numerical results are close to the exact solution, 
y fc w y(tk). The quality of this approximation is measured by the global error, 
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which is defined by ~Eh{tk) = Yk — y(tk)- We say that a method is of order p if 
Ek(i) = 0{hP). 

As mentioned in the introduction, the goal of the analysis presented in this 
paper is to obtain a priori estimates for the global error. The general idea 
is to view the numerical "flow" as a perturbation of the actual flow of the 
equation. As the global error is the difference between the numerical and the 
exact solution, the Alekseev-Grobner lemma can be used to estimate it. 

There is however a problem here: the numerical method provides only values 
at times tk, but not a continuous flow. Iserles and Soderlind (l^] interpolated 
the numerical results to obtain a function defined on the whole time interval. 
Instead, we will use the theory of modified equations, an idea already present in 
the work of Calvo and Hairer || . We note that the same approach is taken by 
Hairer and Lubich jl4|. They proceed to obtain differential equations for the 
various terms in the expansion of the global error. 

The theory of modified equations, analogous to Wilkinson's backward er- 
ror analysis in Numerical Linear Algebra, tries to find a differential equation 
z' = fh(t,z) whose solution is close to the numerical results. Hairer and Lu- 
bich |l3| proved that one can arrange for the solution to be exponentially close 
(in the step size h) to numerical results. This is more than sufficient to apply 
the following theorem. 

Theorem 2. Let y(t) be the solution of the differential equation y' = f(i, y) 
and D<&\ denote its variational flow. Suppose that a numerical method of or- 
der p produces the values {yk}- If the solution z(t) of the modified equation 
z' = f/j(t, z) satisfies z(t/-) — yk = 0{h 2p ) for all k, then 

E h (t) = f D^ s (y(s))d h {s,y^))ds + 0(h 2p ), (5) 
Jto 

where 6 h (t,y) = i h (t,y) - f(t,y). 

Proof. If we take the perturbed equation (|§q) from the Alekseev-Grobner lemma 
to be the modified equation, we deduce 

<t)-y(t)= f D^ s (z(s))S h (s,z(s))ds. 

The left-hand side of this equation is E/i(t) + 0(h 2p ). Furthermore, since the 
method is of order p, we have y(t) — z(t) = 0(h p ). Finally, it follows from the 
theory of modified equations that Sh{t,y) = 0(h p ) (see eg. Jll|). Together this 
proves (||). □ 

Note that the constant implied by the 0(h 2p ) term in (||) depends on t. 
Hence the above result is only valid on bounded time intervals. The same goes 
for the results we will obtain later (Theorems |[ [| and ||). 

To compute the error estimate (|j), we need to know the exact solution y(i), 
the variational flow and the modified equation to compute Sh- Note that 

these can all be determined a priori, hence the above theorem gives an a priori 
estimate for the global error. 
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In fact, if one is interested in the long-term behaviour of the global error, it 
is not essential to know the exact solution over the whole time interval. It often 
suffices to know only the asymptotic behaviour of the solution for large t, because 
this will allow one to compute the dominant contribution to the integral (^). 
This approach is actually taken in the computations that will follow in the next 
sections. 

From now on, we restrict our attention to Runge-Kutta methods. We assume 
that the reader is familiar with the theory of Runge-Kutta methods and B-series, 
as explained in the text books Jl5| and [pOf , amongst others. A B-series is a 
formal power series of the form 

B(a, y) = a(0)y + £ ^ a(r) o(r) P(r) (y). (6) 

Here T denotes the set of rooted trees, a : T — > R is a coefficient function, 
p{r) is the order of the tree r, a(r) denotes the number of monotone labellings 
of r, and F(r)(y) is the elementary differential of the function f associated 
with r. The result of any Runge-Kutta method can be written as a B-series, 
and the coefficient function a depends on the coefficients of the method. 

Hairer |nj uses the concept of B-series to derive a formula for the modi- 
fied equation, taking an approach similar to that of Benettin and Giorgilli 
Suppose that the outcome of the Runge-Kutta method is described by the 
B-series with coefficient function a. Then the modified equation is given by 
z' = i_B(6, z), where the coefficients are recursively defined by 



b{%) = 0, &(.) = 1, and b(r) = o(r) - ]T -^f^r) for p(r) > 2. (7) 

In this formula, db is a Lie derivative, defined as follows. Given two coefficient 
functions b and c with 6(0) = 0, the Lie derivative of B(c, y) with respect to the 
vector field B(b,y) is defined as -^B(c,yb(t)), where yb(t) is a (formal) solution 
of the differential equation y'(t) = B(b,y(t)). It turns out that this is again 
a B-series, and its coefficient function is denoted by dbC. An explicit formula 
for dbC is given in jllj]. 

We now use this expression for the modified equation to rewrite the a priori 
estimate given in Theorem ||. 

Theorem 3. Let y(t) be the solution of the differential equation y' = f(i, y) 
and let D<&\ denote its variational flow. Suppose that this equation is solved by 
a Runge-Kutta method with B-series coefficient function a. Let b be defined as 
in (Q). If the numerical method has order p, then the global error satisfies: 



E h (t)= h^b{r)^-l T {t) + 0{h^) 

where l T {t)= / D$* (y( S ))F(r)(y(s)) ds. 



Here T>2 denotes the set of trees with order at least 2. 
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Proof. As mentioned above, the right-hand side of the modified equation is 
j^B(b,z), with b as in (R). Now F(») = f, so the first term in this B-series is 
the right-hand side of the original equation (|]J). Hence their difference d h (t,y) 
equals 

We now substitute this expression in (§) and move the scalar factors out of the 
integral (remember that the variational flow D& s is a linear operator). This 
results in (|). □ 

Note that the error estimate (||) nicely separates the problem and the method. 
The method only enters the estimate via the coefficients b(r). On the other 
hand, the value of the integrals X T is completely determined by the particular 
equation one is solving. As their role in the global error estimate (||) is similar 
to the role of the elementary differential F(r)(y) in the local error, we will call 
1 T the elementary integral associated with r. 



3 The Airy equation and other linear oscillators 



Iserles studied in detail the leading term of the global error committed by 
various methods, including Runge-Kutta methods, for the Airy equation 



y" + ty = 0, t>0, 2/(0) =2/0, y'(0)=y' . 



0) 



Here we will use Theorem [| to extend his analysis to include higher-order terms. 
This will prepare us for the study of nonlinear variants of the Airy equation, 
viz. ©. ' 

The asymptotic solution of ([)]) for large t can be found by the Liouville- 
Green approximation p^ |, which is also known, particularly to physicists, as 
the WKB(J)-approximation: 



y(t) » A(t)R{9(t))s , where y(t) = 



y(t)' 
y'(t) 



(10) 



Here Sq is a vector whose value is determined by the initial conditions, and 
0(t) — |t 3 / 2 . Furthermore, 



A(t) 



t-v< ' 
t 1 /* 



and R(9) = 



cos tt 
— sin( 



sin 9 
cos 9 



Note that A(t) is a scaling matrix and R{9) is a rotation. 

We want to apply (^) to derive an estimate for the global error, so we have 
to calculate the elementary differentials F(r)(y). This has already been done 
by Iserles ]lS| . He concludes that amongst all the trees of order p, the tall 
tree without branches dominates. We denote this tree by rj (anticipating our 
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analysis of the a nonlinear variant of (f|), compare with the trees in Table 
The elementary differential associated to this tree is F(r p ')(y) = F p y, with 



(-ty o 
o (-ty 



and F 2r +i 



o (-ty 
H) r+1 o 



(11) 



The next step is to compute the elementary integral X \ , defined in (||) . We will 
use the shorthand X p for this elementary integral. Hence we need the variational 
flow D<&\. But since the Airy equation (|])_is linear, the variational flow equals 
the flow map so we can deduce from (|10| ) that 

£>$* = *\ a A(t)ii(0(f) - 0( S ))A- 1 ( S ). (12) 

Now we need to substitute (|l0|), (pd|), and ( |l2] ) in the elementary integral 
X p = J Q D^ t s F p y(s) ds. For even p we find 

k(t)R(9(t) - e(s))A- 1 (s) ■ (-s) r A( S )R(9(s))so ds 

\-s) r ds ■ A(t)R(6(t))s = (-iyf^y( t ). 

r + 1 

For odd p, the calculation is only a bit more complicated. First we compute the 
matrix product that appears in the middle of the integrand: 

A- 1 ( s )F 2r+1 A( s ) = (-l) r s r+ iR(±n). 

Note that R(\n) = o] • With this result, we can tackle the integral, 

\-l) r s r+L 2 ds ■ A(t)R(6(t) + i7r)s = (-l) r ^y R (t). 



r + — 

1 ~ 2 



Here y R (t) is the (asymptotic) solution of the Airy equation (||) with opposite 
phase. It is given by 

y R (t) = A(t)R(6(t) + ±ir)s . (13) 

To find an estimate for the error, we have to substitute the above integrals 
in M). This gives 



E h (i) w dM 2 y(t) + C 2 h 2 t b ' 2 y R (t) + C 3 h 3 t 3 y(t) 

+ C 4 hH 7 / 2 y R (t) + ■■■ + 0(h 2 P), where C k 



(fc + l)!(ifc + l)' 



Note that for a method of order p, the B-series coefficient of any tree whose 
order does not exceed p vanishes. Therefore, the h k term vanishes as well for 
k < p, in agreement with the general theory. 

Note furthermore that the nature of the particular Runge-Kutta method 
used, has no influence on the growth rate of the global error. It does not 
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matter whether the method is for instance A-stable or symplectic. There is one 
exception to this, as pointed out by Orel for methods with b(Tp +1 ) = 0, 
the dominant term in the above estimate vanishes. Hence such methods should 
perform better than general Runge-Kutta methods. 

It is possible to generalize the asymptotic solution (|l^) to a class of nonau- 
tonomous oscillators. Under the assumptions given by Iserles and repro- 
duced in Theorem|| below, the asymptotic solution of the equation y"+g(t)y = 
is still given by (|lf|), except that now 8(t) = f Q \J g(s) ds. The rest of the com- 
putation stays essentially the same. Its result is given by the following theorem. 

Theorem 4. Suppose we are solving the nonautonomous linear oscillator 

y" + g(t)y = 0, t > 0, y(0)=y o , y'(0) = y' o , 

where g(t) is positive for sufficiently large t, \g^(t)\ — o^git) 1 /^ as t — > oo, 
and limsup^oo g(t) = +oo. // we use a Runge-Kutta method with B-series 
coefficient function a, and b is defined as in (0), then the global error is 



E fc (t) 



p<2r<2p 



Er i\r+l b ( T 2r+2) , 2r+l [' 
1 ' (2r + 2)! L 

p<2r+K2p Jt ° 



y R {t) 



g{s) r+ 'ds 



r{t)+0{h 2 V). 



Here p denotes the order of the method, and y(t) and Yji(t) are given by ( [10|) 
and ( |l3| ) respectively, with 8(t) = J* \J g(s) ds. □ 



4 The Emden— Fowler equation 

The previous section should be considered as a mere warm-up to the real chal- 
lenge: nonlinear oscillators. We will consider nonautonomous, nonlinear oscil- 
lators whose behaviour is described by the Emden-Fowler equation 

y" + tV l = 0, t>0, y(0) =yo, 2/(0) = 2/o- (14) 

This equation with v = 1 — n (in which case it is commonly called the Lane- 
Emden equation) was originally proposed to model stars. Considering the star 
as a radially symmetric gaseous polytrope of index n in thermodynamic and 
hydrostatic equilibrium, the relation between the density and the distance to 
the centre satisfies the Lane-Emden equation The Emden-Fowler equation 
also arises in the fields of gas dynamics, fluid mechanics, relativistic mechanics, 
nuclear physics, and the study of chemically reacting systems (see [^7| and 
references therein). 

Note that the choice v = n = 1 reduces the Emden-Fowler equation to the 
Airy equation studied in the previous section. From now on we will assume that 
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n is an odd integer above 1, and that v > — h(n + 3). These conditions assure 
that oscillatory solutions exists, as explained in the survey by Wong We 
remark incidentally that this remains true for any real n > 1, if we replace the 
equation ( |l4| ) by y" + t l 'sgn(y) \y\ n = 0, where sgn(y) denotes the sign of y. 
However, we then lose analyticity at y = 0. 

The equation ( fl4| ) has an obvious scaling symmetry, which can be used to 
reduce the order. Here we will use a different approach though, because the 
asymptotic solution as t — > oo will suffice for our purposes. From now on, we 
set 

v 

n + 3 

The conditions on n and v imply that 7 > — i . Now consider the transformation 
given by 

y{t) = (1 + 2 1 f' { - n - 1 h-"<u{t 1+2 "<). (15) 
If we apply ([l5]) to the differential equation (|l4|), we get 

(1 + 2 1 ) 2n l {n - 1 h^{u"{t 1+ ^) + u n (t 1+2 ^)) 

+ 7(1 + 7 )(1 + 2 7 ) 2 /("- 1 ^- V* 1+27 ) = 0. 

For large t, we can (hopefully) neglect the last term as 7 > — 5, and the above 
equation reduces to u"+u n = 0. The expression + ^(u') 2 is an invariant 

of this equation; it is only an adiabatic invariant to the original equation. We 
will denote the solution of u" + u n — which satisfies the initial conditions 
u(0) = and it'(0) = 1 by w n (t), and note for further reference that this is an 
odd and periodic function. The general solution of u" + u n = is then given 
by u(t) = c 2 /( n+1 ) w n (cit + c 2 ). Note that C\ determines the amplitude of the 
oscillation, while c 2 determines the phase. In other, more sophisticated terms, 
(ci,C2) are action- angle coordinates of the Hamiltonian system corresponding 
to the Emden-Fowler oscillator (see eg. [0). 

It follows that the solution of the Emden-Fowler equation ( |l4| ) is asymptot- 
ically given by 

y(t) « (1 + 27) 2 /("- 1 ) C 2/(n - 1) t-^„( Cl t 1 +^ + c 2 ). (16) 
We repeat that we assumed that n > 1 is an odd integer and that v > — ^(n+3). 

5 Global error of Runge— Kutta methods 

In this section, we repeat the computation of Section || for the nonlinear oscil- 
lator (|l4|). We will assume that the asymptotic solution ( |ilf ) is in fact exact. 
The numerical experiments reported in Section ^ will show that this is a valid 
approximation. 
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The first task is the computation of elementary differentials. For this we 
need to convert the differential equation to a system of autonomous first-order 
equations, 

y[ = 2/2, 

2/2 = -VaV?, (17) 
2/3 - 1- 

Here j/i, 1/2, and j/3 correspond to y, y', and t in the original equation ( p"4| ) 
respectively. 

It follows that the first component of the elementary differential F(r)(y) 
satisfies the following recurrence relations (where the argument y is deleted) 

*!(•)= Ito, F 1 ([)=F 2 (t), F^\y k )=0 (k > 2). 
The second component satisfies 

F2 (V) = -(^f a_kFl(T|) --' fl(rt) (fc - n) ' 

MX/) = ° (fe>n + l), 

where the derivatives with respect to 1/3 are neglected, because they lead to 
terms which are smaller by a factor t 1+21 . Finally, the third component F3 
always vanishes, except for F^*) which equals one. Hence we will drop this 
component from now on. 

It follows from the recursion relations that the first component of the el- 
ementary differential associated with the tree r can be computed as follows. 
Each vertex contributes a factor to F\ (t) . Define the height of a vertex to be 
the distance to the root. Then for vertices whose height is even, the factor is 1/2 
if the vertex is a leaf, 1 if it has one child and zero otherwise, killing the whole 
elementary differential. For vertices with odd height, a vertex with d children 
contributes a factor (apart from a constant) if d < n, and zero otherwise. 

So Fi(t) vanishes if any vertex with even height has more than one child, or 
any vertex with odd height has more than n children. If this is not the case, 
let d' denote the number of vertices with odd height, and p the total number of 
vertices, ie. the order of the tree. Then there are p — d! vertices with even height, 
who have d! children, namely the d! vertices with odd height. Since every vertex 
with even height has at most one child, we have p — 2d' leaves contributing a 
factor 1/2 each and d' vertices with exactly one child, contributing a factor 1. So 
the total contribution of all the d! vertices with even height is i/2~ 2d ■ Similarly, 
the d' vertices with odd height have p — d' — 1 children (the root is not a child 
of any vertex) . So they contribute a factor t d v yl A ^ d . We conclude that 

fi(r)(y(t)) - C 1 , r 1? v vt*' 1) *~'* 1 (t)v$-* e (t). 
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By a similar reasoning, wc find that 



F 2 (r)(y(t)) = C 2 , T &- d >y^- {n+1)d '(t)y 2 2 d '- p+1 (t). 



Note that the conditions on the tree which make sure that the elementary dif- 
ferential does not vanish, imply that all the exponents in the above formulae 
are positive. Also, the only trees for which both F\ and F% are nonzero, are the 
tall branchless trees r ' . 

If we substitute the approximate solution (^) , we find that the elementary 
differentials are 



F(r)(yW) 



' C 3 , T t^ 2 P-Vw { r? +1)d '- p+1 (i) w' n p - 2d ' (t) ' 

c 4iT t^ 2 p+Vw7- (n+1)d \t)w' n 2d '- p+1 (t) 



(18) 



where t = cit 1+ ' 21 + c 2 . The growth rate of the elementary differential is de- 
termined by the exponent of t. Note that the variable d! does not enter in this 
exponent. The surprising conclusion is that all trees with the same order con- 
tribute a term with the same growth rate, independent of their shape. This is in 
stark contrast to the linear Airy equation, where the differential corresponding 
to the tall tree r p l dominates. 

The next step is to calculate the elementary integral X T . For this, we need 
to multiply the above differential with the variational flow matrix and integrate 
the resulting expression. To compute the variational flow, we introduce the 
map Xt : R 2 — > R 2 , defined by 



(1 + 2 7 ) 2 /("- 1 ) C 2/( "- 1) t-^„(c 1 i 1 +^ + c 2 ) 
(1 + 2 7 ) 1 + 2 /("- 1 ) c ; +2/( "- 1 V<(c 1 t 1 + 2 ^ + pa) 



(19) 



So X t maps the parameter space to the solution space at time t. Neglecting 
lower-order terms, it follows that the flow map satisfies <E>* = X t o X~ l . Hence 
we can write the elementary integral (|^) as: 

Ir(t) « DX t (c) [ CA s - 1 (y( S ))F(r)(y( S ))d S . (20) 
Jo 

To find the integrand in the above expression, we multiply the inverse of the 
Jacobian matrix of ( |l9| ) with (18). The result is 

s 2 ^[c^w^ +1){d ' +1) - p {w' n y- 2d ' + c & , T w7- (n+1)d '«) 



/ \2d' -p+2 



„27p+2'y+l ( n r,.,( n + 1 )(- d ' + 1 )-P(r,.,i \P-2d' , r< ,„ np 
5 U ' Iw,tW„ ( w n) + ^8,rW n 



-( n+1 ) d ' f w > )2c2'-p+2 



(21) 



where the functions w n and w' n are evaluated at s = cis 1+27 
In the calculation, we used that 



c 2 . 



<(t) = -<(«) and < 2 (i) = - H i T <+ 1 (f) + l. 
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The first equality is indeed the definition of w n , and the second one follows 
by multiplying the first one by w' n (t) and integrating, using the initial condi- 
tions io n (0) = and w' n (0) = 1. 

The next step is to integrate (|2l|). But consider the exponents of w and w' . If 
p is even, these exponents are also even and hence the integrand is nonnegative. 
Now recall that w is odd and periodic. So in the case when p is odd, the 
functions w and w' are raised to an odd power, which means that the integrand 
oscillates around zero. Thus we can expect cancellations in the latter case, but 
not if p is even. We stress that this phenomenon does not occur in the linear 
case, analysed in Section [| 

In fact, we have 



wi(s)w' n m (s) ds 



if either i or to is odd, 
if both £ and to are even. 



Here Ce mn (t) denotes a oscillatory function with the same period as w n {t), and 
Czmn is a constant. It follows that 



s k w l n {s) w' n "\s) ds 



c kemn (t)t k -^ + o(t k -^- 1 ) 
C Mmn t k+1 +o{t k -^), 



if i or to are odd, 
if £ and to are even, 



where again Cktmn and Cktmn denote a periodic function and a constant, re- 
spectively. 

We can use this result to integrate (|2l]), and find 



f DX- 1 F(r)(y)ds= < 
Jo 



C 9 . T {i) t 2 1P- 2 1 + 0(*27P-47-l) 

c 9 , T t 2w+1 + o(t 2 ~<p- 2 ~<) 

Cio,^ 27p+27+2 + 0{t 2 ~'P +1 ) 



, if p odd, 



if p even. 



(22) 



To compute the elementary integral X T , we need to premultiply the integral (|22|) 
with DX t , cf. (|2^). But the expression (^2|) has an interpretation by itself. 
Remember that X t maps the parameter space to the solution space. So the 
integral ( |22| ) represents the error in the parameter space. We conclude that 
the energy error associated with the tree r grows as t 2lp+1 if p is even, and 
as t 2l9 ~ 21 if p is odd. The second component gives the phase error. 

Multiplying the Jacobian matrix of the map X t with the integral (|2^) gives 
us the elementary integrals 



l T (t) 



Cii, T {t)t 2 iP-~<+ 1 + 0{t 2 iP- i ~t) 

C 12 ,r(t) t 2 ~<P+~< +1 + 0{t 2 ~*P-~<) 

' Cii, T {t)t 2 -<p+- <+2 + o{t 2 ~<p-~< +l ) 

Cl2, T [t) t 2 "IP+ Z ~<+ 2 + 0(t 2 W+7+l) 



if p odd, 



if p even. 
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Finally we can find an estimate for the global error by adding the contributions 
of all trees according to (0). For the first component, we find 

E h (t) m C 1 {t)ht^ +2 +C 2 (Wt^ +1 +C 3 (t)hh 9 ^ +2 +C i (t)hH 9 '< +1 + - • -+0(h 2p ). 

More formally, we have the following result. 

Theorem 5. Suppose one is solving the Emden-Fowler equation y" + t u y n = 0. 
with v > — i(n+3) and n an odd integer greater than 1, with a numerical method 
of order p that can be expressed as a B-series. Then the global error has the 
form 



E h (t)= ]T h 2r 

p<2r<2p 

+ h2r+1 

p<2r+K2p 

Here C\(t) denotes a function periodic ini — ci£ 4 / 3 + C2, and 7 = v/(n + 3). □ 

We would again like to draw the reader's attention to the difference between 
the even and the odd powers of h. Unfortunately, we do not know what causes 
this phenomenon, and whether it also occurs for other differential equations. 



r 1 

L/ 2r 



(t) £ 4 7>-+7+l + 0^4 7 r- 7 ) 

CL(t) t^+ST+i + o(t^ T +"f) 



C2r+lW t^ r +^'+ 2 + O (£47^+37+1) 
Cfr+lft) t^ r+7 ^+ 2 + 0(i 4 7'-+57+l) 



0(h 2p ). 



6 Numerical experiments 

The purpose of this section is to supplement the calculations of the previous 
section with some numerical experiments. In particular, we want to see whether 
we were justified in using the asymptotic solution (^). Furthermore, it could be 
interesting to study whether the 0(h p ) term of the global error E^(t) dominates, 
or whether other terms have to be taken into account. 

All experiments are performed with the parameters n = 3 and v = 1, so the 
differential equation we are solving is 

y" + ty 3 = 0. (23) 

The reason for this particular choice is that the function w n (t), the solution of 
u" + u n = satisfying u(0) = and u'(0) = 1 (cf. Section |), can be expressed in 
terms of Jacobi elliptic functions (see eg. In fact, we have w 3 (t) = sd(t \ i). 

The parameter i will be dropped from now on. As a consequence, we can 
calculate the elementary integral associated with any given tree explicitly. For 
the first couple of trees, this yields the results listed in Table [l]. 

Our first example concerns Runge's second order method, given by 

y«+i = yVi + hf(t n + ~h,y n + lhf(t n ,y n )). (24) 



13 



Tree 
r 


Elementary differential Elementary integral 

F(r)(y) J T (t)=/ t £>$ t s F(r)(y)d S 


A: I 




-ViVa 

r^vlv^ ~ vl. 




0(i 5 / 6 )_ 




128 AT 4 + 17/6 -j'/TS" 

_|§V2cf x t^/Sd 3 (t)_ 










-63/12/22/3 - §y\v2 






0(i 7 / 6 )_ 




"-^^cf X * U/6 sd'(t)" 
_||V2cfxt 13/6 sd 3 (f)_ 




4 I 








'C(< 5 / 6 )" 

0(< 7 /6)_ 




-fiV2c! X i 11/6 sd'(f)" 
_i|V2cf X ^/ 6 S d 3 (f)_ 






[ 









0(i 3 / 2 ) 




J^V2c 7 i 23 / 6 sd 3 (f)_ 




*\} 






pvtyivl + 3y 1 ?/3 






0(i 3 / 2 )_ 




-^>/2 sd 8 (t) 




rt-Y 




-^vwlvi - 6yiv2 






0(i 7 / 6 )~ 





"-Sv / 2^ 7 / 2 sd'(f ) " 








9yjWI + 3yi2/3 




"0(t 7 / 6 )~ 
0(i 3 / 2 )_ 




liiV2 C ?^ 2 sd'(f) 
-^V2c^/Sd 3 (t) 





Table 1: Trees of order < 4, with their elementary differentials and integrals for 
the Emden-Fowler equation (p3|). The third component is suppressed as it is 
always zero. In the last column, x — jk Jo s d 2 (s) ds where 4K is the period of 
the function sd, t = cii 4 / 3 + C2, and only the term of leading order is displayed. 
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This is a Runge-Kutta method, hence it can be written as a B-series B(an,y). 
Some values of the function an are 



ok(«) = i» «k(I) = i> ««(V) = i «fl(NI/) : 

=« K (Y)=o- 



It follows from (^) that the B-series coefficients of the modified equation are 
&«(•) = !. MI) = 0, 6k(V)=-j. Mi)=" 1 . 



6ji(NI/)=0, 6 fl («Jj=0, bj- 



= 3, 



,(Y 



If we use these numbers in the error estimate (JSJ) , together with the elementary 
integrals from Table [l], we obtain the following estimate for the global error 
committed by Runge's method: 



^V2cf^ 2 sd'(f) 

512 



1871 



T V2clt 23 / 6 sd 3 (t) 



(25) 



To check this estimate, we compute the solution of the nonlinear oscillator ( |23j ) 
with Runge's second-order method (|2J). The initial conditions are y(0) = 1 and 
y'(0) = 0, which lead to a solution with c\ « 0.7. The solution is compared 
to the result of the standard fourth-order Runge-Kutta method with step size 
h = 1/10000. According to Theorem ||, this would give an error of about 10~ 9 , 
so we can consider this to be the exact solution. Hence the global error Eh(t) 
is computed by subtracting the result of Runge's method from the "exact" 
solution. The first component is depicted in Figure [i] The left column shows 
the time interval [0, 50], and on the right the larger interval [0, 2000] is displayed. 

In the left column of Figure [l], the oscillating curve shows the global er- 
ror Eh(t). The oscillating behaviour is correctly predicted by the error esti- 
mate (p5|); this is the factor sd'(t). The amplitude of the oscillations, as pre- 
dicted by ( p5j ) , is shown by the thick curve in the left-hand column in Figure |l|. 
We see that the estimate (25) describes the actual error accurately. 

The right-hand column of Figure [l] shows a much larger time interval. Here 
the oscillations are compressed so heavily that the error appears as a grey blob. 
The dashed curve shows the first, leading term of the error estimate (p5|), and 
the solid curve shows the sum of both terms. We conclude that the leading 
h 2 term of the estimate does not describe the actual error correctly, but that 
the error is predicted accurately if the h 3 term is included. For h — 1/1000 
the latter estimate breaks down around t = 1200. We note that at that point, 
the amplitude of the solution is approximately 0.3 and the error has about the 
same size, so the numerical solution has deviated considerably from the exact 
solution. For smaller values of the step size, the estimate ( p5| ) is accurate over 
the entire time interval [0, 2000]. 
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Figure 1: The first component of the global error committed by Runge's second- 
order method, and the estimate (p5j). 



1G 



It follows that for large t (which here means: in the order of 1000), Runge's 
method behaves essentially as a third-order method. To check this, we compare 
it with Heun's classical third-order method: 



ki = f(t„,y„) 
k 2 = f(t„ + |/i,y„ 
k 3 = f(t n + §/i,y„ 
Yn+i = Yn + h(l\s.i + 



Butcher tableau: 



|k 3 ) 








1/3 


1/3 


2/3 


2/3 




1/4 3/4 



(26) 



A similar calculation as for Runge's second-order method gives the following 
estimate for the global error committed by this method: 



E h (t) 



512 
3572 
1024 
107163 



T V2 4t 7 / 2 Sd'(i)' 



50208 



229635-V2cfi 5 / 2 sd'(i) 
60416 



V2clt 17 / e sd 3 (i) 



688905 

557056 /o „8 £25/6 j> (I\ 
34543665 VZC 1A C hu \ L ) 

1114112 . AT„9.,i9/2„j3/7- 
103630995 



V2c? X t 9/2 sd 3 (t) 



(27) 



The actual error and the above estimates, for step size h — 1/2000, are dis- 
played in the second row of Figure |^. The first row displays Runge's second- 
order method. We see that the error estimate ( p7| ) again provides an excellent 
description of the actual error. Furthermore, the difference in order between 
Runge's and Heun's methods shows clearly for small values of t (see the left- 
hand column in Figure ^). For large values of t however (cf. the right-hand 
column), Runge's method behaves essentially as a third-order method and we 
see indeed that the difference between the two methods is much smaller. 

The bottom row in Figure [2] shows a specially tuned method. Remember 
from Section [| that for the Airy equation (^) , the contribution of the tree Tj] to 
the global error dominates, and this is used by Orel [£3] to construct specially 
tuned Runge-Kutta methods. We found in Section || that this is not the case for 
the Emden-Fowler oscillator (14). However, if we carefully study the elementary 
integrals T T in Table |], we see that they are scalar multiples of each other: 

4X T «(t) = -12J Tl (t) = \21 Tl {t) = -3T j(t). 

Therefore, bearing in mind the factor a(r) in Theorem ^, the h 3 term will be 
killed for a third-order method with 



36(\lV)-3&(^) + fc(Y)-4& 



0. 



From the general theory, it follows that a method has order three if 
&(•) = ! and 6(1) = b(\/) =&($) = 0. 

These are five conditions together. An explicit 3-stage Runge-Kutta method 
has six free parameters, so there is some hope that we can find such a method 
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Figure 2: The first component of the global error committed by Runge's second- 
order method fl24|), Heun's third-order method (p6|), and the specially tuned 
third-order method (^8|), all with step size h = 1/2000, together with their 
respective error estimates (25|), (|2~7|), and (29). The lines in the plot have the 
same meaning as in Figure 



18 



satisfying these five conditions. Indeed, it turns out that there exists a one- 
parameter family of these methods. A particular instance is: 



ki 
k 2 



f(t n + h,y n 



Butcher tableau: 



hki) 



k 3 = f(t„ + §/i,y„ + ffcki - |/ik 2 ) 



yn+i — y« 










1 


1 




3/2 


9/4 


-3/4 




7/18 


5/6 -2/9 



(28) 



If we compute the estimate for the global error of this method, we find 



E„(i) « /i 4 



5008 



_ >: , 51 ,^ C ?tV 2sd '(i) 

iiiV2c^ 17 /6 sd 3 (f) ^ 



78848 



1279395 
157696 
3838185 



y2 c 8 xt 25/6 gd '(£y 

y2c? X t 9/2 sd 3 (f) 



(29) 



As expected, the ft 3 term has disappeared. This estimate, together with the 
actual error, is displayed in the bottom row of Figure ||. We see that the 
global error of this method is much smaller than that of Heun's method, though 
both methods are of third order. The other thing to note is that the error 
estimate (^9|) is not as accurate as for the other methods, even though it predicts 
the right order of magnitude. 

We stress that we are not recommending the use of the method (p8h ; it works 
well for this equation but we know of no reason whatsoever why it should also 
perform well for other equations. It may however be important to gain a better 
understanding of the mechanisms that cause this method to work so well. 

The reader should keep in mind that the local error of the method 
is still 0(h 4 ), as for all third-order methods. In fact, the B-series coefficient 
function ax of the method (E8) satisfies 



So the local error of this method is 



B(a T - e,y) 



Ij/iWI + aJ/iJ/3 + Ivhz + \yivl 



o. 



0{h b ). 



where e denotes the function with e(0) = and e(r) = 1 for r ^ (the B-series 
B(e,y) gives the exact solution (II), so B(a,T — e,y) is the local error). 

It follows that the global error Efc(t) must be 0(h 3 ), in contrast to the error 
estimate (|9|). It is only the term of order h 3 [t 7 ^ 2 , t 23 ^ 6 ] T , which for generic 
third-order methods dominates the h 3 term, that disappears. But there still 
is an h 3 contribution to the global error. We believe that this contribution 
explains the substantial difference between the estimated and the actual error 
for this method, as displayed in Figure |[ The difference in the shape of the 
oscillations of the global error for the last method supports this conclusion. 
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7 Discussion 



We used a combination of the Alekseev-Grobner lemma, the theory of mod- 
ified equations, and asymptotics to calculate an estimate for the global error 
committed by a Runge-Kutta method for a class of linear and one of nonlin- 
ear oscillators. Numerical experiments show that these estimates are generally 
accurate. 

While in the linear case, the global error behaves more or less as we expected 
after the analysis by Iserles , the nonlinear case brings some surprises which 
indicate that the entire calculation of Section || is necessary. The specially 
tuned method ( p8| ) shows that it is not enough to consider only the local error. 
It neither suffices to study just the first term of order h p ; Figure || shows that 
the next term sometimes dominates. Finally, a comparison of Theorems |] and |^ 
shows that linear and nonlinear oscillators behave differently. 

Let us highlight these differences. In the linear case, the contribution of one 
single tree (namely the tall, branchless tree) dominates the local error. This 
is no longer the case for the nonlinear oscillator (|l4|). On the other hand, the 
contributions of all trees of the same order to the global error are constant 
multiples of one another (this has been checked for order < 6). This allows us 
to kill the leading term of the global error. Experiments show that this strategy 
indeed works, but we still do not completely understand what is going on. 

The other difference concerns the integral in the error estimate. The integra- 
tion may bring down the time exponent due to oscillations in the integrand. For 
our equation this happens for even powers of h. As a consequence, the lowest 
order term in the expansion of the global error in powers of h does not always 
reveal the true behaviour of the global error, as the next term may dominate 
for sufficiently large values of t. 

We stress that we have only considered one particular class of nonlinear 
oscillators, namely those satisfying the Emden-Fowler equation (|l4|). Hence we 
do not know how general these phenomena arc. We plan to try and repeat the 
above analysis for other equations. An obvious target is the generalized Emdcn- 
Fowler equation, given by y" + a(t)y n — with a(t) > (27j. A more ambitious 
plan is to try and connect our analysis to the recent interest in highly-oscillatory 
systems, triggered by Garcia- Archilla, Sanz-Serna and Skeel ||, and Hochbruck 
and Lubich 0. 

A delicate issue is the validity region of the estimates derived in this paper. If 
t is too large the 0(h 2p ) term will probably dominate, rendering the estimates 
worthless. On the other hand, the estimates of Sections ||, ||, and ^ use the 
asymptotic solution as t — > oo, so we need t to be sufficiently large. In any case, 
the step size h needs to be sufficiently small; but if it is very small only the 
hP term will contribute. The final thing to keep in mind is that the expansion 
of the modified equation in powers of h usually diverges. The only definitive 
statement we can make is that more research needs to be done on this issue. 

This paper has concentrated on Runge-Kutta methods. Of course, there 
are many other methods. If one knows that the solution will be oscillatory, one 
is probably better off with a method that is specifically designed for this case. 
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The reader is referred to the review paper by Petzold, Jay and Yen Q for 
examples. We hope that the framework of Section |^ is sufficiently flexible to 
include a wide variety of methods. This requires a theory of modified equations 
for these methods, as Hairer and Faltinsen H] have done for multistep and 
Lie-group methods respectively. Another challenge is presented by variable step 
size methods. 

Let us end by comparing the general estimates from Section |L especially 
Theorem 2|, with some results from the literature. The estimate derived by 
Iserles |L8| can be written as 

Eft(t) = i jf D^ s (y(s)) L h (s, y(s)) ds + 0(^ +1 ), (30) 

where Lh(t,y) denotes the local error of the method. Indeed, this estimate 
follows from (||). We note that ( ^0|) gives only the 0{h p ) term of the global 
error, and we know that this does not always paint an accurate picture. On the 
other hand, the local error is easier to determine than the modified equation. 

Earlier, Viswanath ]2(| took quite a different approach but derived a bound 
on the global error which can be brought to almost the same form as (|30|), 
namely 

E h (f) < / ||^*(y(s))|| ds ■ ~ max \L h (t, y(t))\ + 0{W +1 ). 

This bound severely overestimates the actual error in some circumstances. On 
the other hand, Viswanath [M showed that the first factor can be estimated 
for some large classes of problems, yielding an upper bound for the global error 
of any reasonable method. 

The final note should be on the older but still highly relevant work of Hairer 
and Lubich who built upon Gragg's asymptotic expansion ||. They ex- 
panded the global error in powers of the step size h, and showed how the terms 
can be found recursively by solving a differential equation involving the previous 
term. In principle, all terms of the asymptotic expansion of the global error can 
be found in this way. Unfortunately, we are not able to solve the differential 
equations that appear when we apply their approach to the nonlinear oscilla- 
tor (^). The first p terms of Hairer and Lubich's expansion must obviously 
equal the integral in (^). However, the connection has not been found yet. 
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